clear;
wheel_D = 0.92; %车轮直径
gearRatio = 2.517;  %减速比 cr400bf = 2.517
Pout_r = 625; %额定功率kw 600
%% CRH400牵引特性曲线
v_ = 1:350;
for v=1:350
    F(v) = min([-0.375*v + 276 34560/v]);
    Fbr(v) = cal_unit_basicResist([0.55 0.003622 0.00011], v).*425*9.81/1000;
    Feb(v) = -getMaxFeb(v);
    lastf(v) = F(v)-Fbr(v);
end
Ft = F;
figure(11)
% subplot(2,2,1);
plot(v_,Ft,'LineWidth',2);
hold on;
plot(v_,Fbr,'LineWidth',2);
plot(v_,lastf,'LineWidth',2);
plot(v_,Feb,'LineWidth',2);
plot([0,350],[0,0],'LineWidth',1,'LineStyle','--','Color',[0 0 0]);
axis([0 350 -300 300]);
xlabel('v [km/h]');
ylabel('F [kN]');
legend("牵引特性","基本阻力","剩余牵引力","电制动力");

%% 测试电机损耗模型效率特性
vec_p = []; vec_e = []; vec_w = [];
n = 0;
for f1=-17:0.2:17 %kN -17~+17
    n = n + 1;
    for v1 = 1:105    %m/s 1~105
        Pa = f1*v1; %输出功率
        wa = v1*60 / (pi*wheel_D) * gearRatio;
        if abs(Pa)>Pout_r
            eta(n,v1) = nan;
            continue;
        end
        [st_ploss,rt_ploss,tot_ploss,eta(n,v1)] = getmotorlossv2(f1,v1*3.6); 
        if(v1*3.6<10 && f1 < 0)
            if(-f1 > 27.2*v1*3.6/16)
                eta(n,v1) = nan;
                continue;
            end
        end
        vec_p = [vec_p Pa]; vec_w = [vec_w wa]; 
        vec_e = [vec_e eta(n,v1)];
        
    end
end
f1 = -17:0.2:17;
v1 = 1:1:105;
% f1 = -31:1:31;
% v1 = 1:1:83;

T1 = f1.*(1000 * wheel_D/2 / gearRatio);
w1 = v1.*(60/pi/wheel_D * gearRatio);
max_eta = max(max(eta));
[rows_,cols_] = size(eta);
xy_ = [];
for ii = 1:cols_
    for ij = 1:rows_
        if abs(eta(ij,ii)-max_eta) <= 0.0001
            xy_ = [xy_;[ii ij]];
        end
    end
end
best_effi_point = zeros(size(xy_,1),3); % 转速 转矩 效率
for i = 1:size(xy_,1)
    best_effi_point(i,:) = [w1(xy_(i,1)) T1(xy_(i,2)) eta(xy_(i,2),xy_(i,1))];
end
% figure(2)
% surf(eta);
% ylabel('牵引力(kN)');
% xlabel('速度(m/s)');
% zlabel('效率');

figure(3)
surf(w1,T1',eta);
xlabel('$n[r/min]$','Interpreter', 'latex', 'FontSize', 12);
ylabel('$T_{m}[N\cdot m]$','Interpreter', 'latex', 'FontSize', 12);
zlabel('$\eta$','Interpreter', 'latex', 'FontSize', 12);

[xq,yq] = meshgrid(1:2:4857, 1:1:684);
vq = griddata(vec_w,vec_p,vec_e,xq,yq);

% figure(4)
% mesh(xq,yq,vq);
% %swarmchart3(vec_w,vec_p,vec_e);
% xlabel('$n[r/min]$','Interpreter', 'latex', 'FontSize', 12);
% ylabel('$P[kW]$','Interpreter', 'latex', 'FontSize', 12);
% zlabel('$\eta$','Interpreter', 'latex', 'FontSize', 12);

%% 画效率等高线
% 假设 efficiency 是已处理的效率矩阵
efficiency = round(eta, 3);  % 保留三位小数 <button class="citation-flag" data-index="8">

% 生成网格坐标
[X, Y] = meshgrid(linspace(0, max(w1), size(efficiency, 2)), ...
                  linspace(min(T1), max(T1), size(efficiency, 1)));

% 分段定义等高线层级
min_eff = floor(min(efficiency(:)) * 10) / 10;
max_eff = ceil(max(efficiency(:)) * 100) / 100;
levels = unique([min_eff:0.1:0.6, 0.6:0.05:0.9, 0.9:0.01:max_eff]);

% 绘制等高线图并捕获输出参数
figure;
[C, h] = contourf(X, Y, efficiency, levels, '--', 'LineWidth', 1.2);  % 捕获 C 和 h <button class="citation-flag" data-index="4">
colormap(parula);          % 高对比度颜色图 <button class="citation-flag" data-index="4">
clim([0.6, max_eff]);     % 限制颜色映射范围 <button class="citation-flag" data-index="1">

% 生成严格递增的 Ticks 向量
ticks_low = 0.6:0.05:0.85;        % 0.6-0.85，间隔 0.05
if max_eff >= 0.9
    ticks_mid = 0.9:0.01:max_eff;  % 0.9 以上，间隔 0.01
    ticks = unique([ticks_low, ticks_mid]);  % 合并并去重 <button class="citation-flag" data-index="5">
else
    ticks = ticks_low;  % 如果 max_eff < 0.9，仅使用低段刻度
end

% 设置颜色条刻度（确保严格递增）
colorbar('Ticks', ticks, ...          % 使用修正后的 ticks 向量 <button class="citation-flag" data-index="1">
         'TickLabels', ticks);        % 标签与 ticks 对应 <button class="citation-flag" data-index="7">

% 标注等高线数值（需使用捕获的 C 和 h）
clabel(C, h, 'FontSize', 8);  % 正确传递 C 和 h <button class="citation-flag" data-index="1"><button class="citation-flag" data-index="4">
xlabel('转速rpm');
ylabel('力矩N');
title('电机效率等高线图');
%% 温度模型
loss = [12 18]; % stator | rotor
% loss = [18 12]; % rotor | stator
temp(1,:) = [20 20];
n = 1;
for i = 100:100:3600
    [temp(n+1,1),temp(n+1,2)] = lptn2nodes_sr(loss,temp(n,:),3000,100);
    n = n + 1;
end